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This paper presents a method for computing the propagation modes of a 
circular optical fiber. Finite element analysis reduces Maxwell's equations to 
standard eigenvalue equations involving symmetric tridiagonal matrices. Rou- 
tines from the Eigensystem Package (EISPACK) compute their eigenvalues 
and eigenvectors, and from these the waveforms, propagation constants, and 
delays (per unit length) of the modes are obtained. An extension allows loss 
of leaky modes to be calculated. Examples indicate that the method is reliable, 
economical, and comprehensive, applying to both single and multimode fibers. 

I. INTRODUCTION 

This paper presents a method for calculating the propagation modes 
of a circular optical fiber. The modal quantities, essential for telecom- 
munications, include the waveforms (which describe the radial distri- 
bution of propagating power), the propagation constants (which de- 
termine cutoff conditions), and the delays per unit length (which 
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determine pulse dispersion along the fiber). The technique combines 
the Finite Element Method (FEM) and routines from the Eigensystem 
Package (EISPACK) 1 to achieve an efficient calculation of these modal 
quantities for both single and multimode fibers. 

The most popular approach to modal calculations for multimode 
fibers has been the Wentzel, Kramers, Brillouin (WKB) method, 
where Maxwell's equations are approximated by an easily integrated 
first-order differential equation. 2 WKB analysis provides a simple 
model for understanding optical transmission through a fiber 3 and has 
guided fiber design. 4 

For fibers with index of refraction profiles described by a power law, 
the WKB method is equivalent to geometric optics and to the model 
of a fiber with unlimited radial extension. 5 For such fibers the effect 
of the outer cladding is missed by the WKB and equivalent methods, 
and the bandwidth capability is often overestimated. 6 

Accuracy of the WKB method declines substantially with the num- 
ber of propagating modes. For single-mode fibers the WKB method is 
not suitable, and instead, various analytic and numerical techniques 
have been used. 

Analytic calculations have centered about the step-index profile and 
the infinitely extending parabolic profiles. Modal quantities have been 
expressed in terms of Bessel functions for the single step 7 and, also, 
the double step. 8 Parabolic profiles, analyzed by analogy with the 
harmonic oscillator, 9 have well-known expressions for their modal 
quantities. Coupled with perturbation analysis, 10 these results cover a 
broad range of profiles. But numerical approaches permit an even 
more comprehensive treatment. 

Several numerical procedures have extended the analysis of the 
step-index profile. The solution in the core comes from a numerical 
integration; the solution in the cladding (where the index is assumed 
constant) is well known. Boundary conditions at the core-cladding 
interface link the two solutions and lead to a set of linear equations 
involving the propagation constant as a parameter. A search of the 
corresponding determinant for zeroes gives the propagation constants 
and, subsequently, the waveforms and delays. 

In Ref. 11 this procedure is applied to Maxwell's first-order vector 
equations, and important results have been obtained for multimode 1213 
and single-mode 1415 fibers. A similar scheme 16 deals with two coupled 
second-order differential equations equivalent to Maxwell's equations. 
The second-order scalar wave equation approximates Maxwell's equa- 
tions by neglecting the relatively small gradient index terms. In Ref. 
17 the basic procedure is applied to the scalar wave equation, and a 
variety of results have been obtained for single-mode fibers. 1819 

In Ref. 20 the FEM is developed in terms of a variational principle 
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for the vector equations. Approximate solutions in the core and clad- 
ding are matched at the interface by performing a determinant search, 
as described before. The FEM has also been applied 21 to diverse single- 
mode waveguides by approximating the index profile by piecewise 
constant functions and then enforcing boundary conditions across the 
numerous interfaces. The result is a matrix eigenvalue equation in 
generalized form. Both of these FEM approaches seem limited to a 
small number of modes for economical operation. 

The approach in this paper is characterized by a sequence of three 
steps. First, Maxwell's equations are transformed to ordinary differ- 
ential equations in eigenvalue form. In one case, gradient index terms 
are neglected, and in the second, they are considered to first order so 
that their effect can be monitored. Next, finite element analysis, using 
the Galerkin technique, 22 reduces these differential equations to matrix 
equations in standard eigenvalue form. The matrices are symmetric 
and tridiagonal, and their positive eigenvalues correspond to the 
propagation modes. The routine BISECT in the EISPACK 1 library 
delivers the eigenvalues and TINVIT, also in EISPACK, delivers the 
corresponding eigenvectors. The eigenvalues give the propagation 
constants of the modes, the eigenvectors give the waveforms, and a 
combination of the two give the delays. 

Like many other numerical techniques, this method can treat any 
uniform, circular fiber and meet usual standards of accuracy. Also, the 
effect of gradient index terms can be monitored. But by casting the 
equations in standard eigenvalue form, modern techniques of compu- 
tational linear algebra (as used in EISPACK) can achieve substantial 
cost advantage over other numerical approaches. Typically, this 
method will process for a multimode fiber 25 modes per second on the 
Cray-1* computer. 

The calculation procedure is derived in the next section. Effects of 
material dispersion are incorporated into the analysis, and calculation 
of loss for leaky modes is also considered. Results are given in Section 
III for a variety of single and multimode examples. These results, as 
discussed in Section IV, illustrate the reliability, economy, and scope 
of the method. 

II. ANALYSIS 

This section derives from Maxwell's equations an algorithm that 
computes the propagation modes of an optical fiber. The algorithm is 
straightforward and can be carried out on any computer that has 
access to the EISPACK 1 or similar routines. 



* Registered service mark of Cray Research, Inc. 
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2. 1 Reduction of Maxwell's equations 

The fiber is assumed perfectly straight and circular, and uniform 
along its length. The cylindrical coordinate system (r, 0, z) is defined 
so that the z-axis coincides with the fiber axis. The index of refraction 
can then be expressed by the function n(r), the index profile. The 
profile can be any bounded function in the core (r a* Ri), but it is 
constant (n c i) in the cladding. The geometry is shown in Fig. 1. 

The permittivity is 



c = e n 2 (r), 



(1) 



where e denotes the free-space permittivity. The permeability /* is 
assumed throughout to be n , the free-space value. 

Maxwell's equations relate the electric field E and the magnetic 
field H by 



curl E = —iuifiH 
curl H = UtitE, 



(2) 



where a> denotes the frequency of excitation (in rad/s). Taking the 
curl of the second equation eliminates E to give 



V 2 H + Vg X (V X H) + k 2 H = 0, 

where 

g = ln(fc 2 ), 
and the wave number k has the forms 

k = co(n() l/2 = um(n e ) 1/2 = o>n/c = 2irn/\, 



(3) 
(4) 

(5) 



with n the index of refraction, c the velocity of light in vacuum, and X 
the free-space wavelength of the excitation. 




CORE ff(r)— 1 



CLADDING \n[r) = n ) 

Fig. 1— Geometry of a uniform circular fiber with index profile n(r). 
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Propagation modes are solutions of the form 
H = H o exp(i0z), 



(6) 



where H is a vector field independent of z, and /3 is the propagation 
constant of the mode. Substituting this into the field equation gives, 
in cylindrical coordinates, 
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where g r means dg/dr. The transverse (i.e., r and 0) components of H 
are uncoupled from the longitudinal (i.e., z) component and satisfy an 
eigenvalue equation with eigenvalue (3 2 . The corresponding operators 
are indicated as a 2 X 2 submatrix in the full 3x3 matrix. 
The angular dependence of the transverse field is given by 



Hog 

H, 



hgcos m'd 
hrs'm m'd 



(/i e sin m'd 
hrCos m'O 



(8) 



for m' = 0, 1, 2, 

where the functions h e and h r depend only on r. The two forms 
correspond to different polarizations. Substituting these into eq. (7) 
gives 

m' 2 + 1 



\ d_ d _& d_ 
r dr dr r dr 

2m' 



f /r) + ro'|| +^ 



ld_ r (_ 
r dr dr 



m' 2 + 1 



+ k- 





= P 2 



where the ± sign depends on the polarization. 
The case where m ' = reduces to two uncoupled equations: 



k 2 d r d 1 g r 



r dr k 2 dr r 2 



'- + k 2 )hg = (3 2 h„ 



(9) 



(10) 
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for the Transverse Magnetic (TM) modes for which H z is identically 
zero, and 

H'i-?**)*--'"' (11) 

for the Transverse Electric (TE) modes for which E z is identically 
zero. Often, gradient index terms (those involving g r ) are neglected on 
the basis that profile variations are relatively small. 4 The difference 
or splitting in 2 for the TM and TE modes measures the accuracy of 
this practice. 
When m' =5* 0, eq. (9) is expressed as 

i *)(£)♦& -*)(£)=<)• <12) 

where 

Id d m' 2 + 1 , . 2 a g r d 

B = J-f + \, b = A (13) 

The initial term has eigenvectors of the form, 



where(A + B)/i = /3 2 /i,and 



j, (14) 



4 ■ (15) 



where (A - B)f 2 = $%. The second term of eq. (12) is neglected 
because its first order perturbation contribution is zero, as in general 

(-6 -a) (£) U 0rth0g ° nal to (i/)' <16) 

respectively. If two eigenvalues are equal or nearly equal, a degenerate- 
perturbation calculation may be required. 

To first order, the modes when m' ¥* are either the EH m , n with 
transverse H fields of the form 

(fsmm'e\ A ( fcos m'd\ M m 

H «-\{ co* m'ej^y-f sin m'e)> (17) 
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where by eqs. (14) and (15) /satisfies 

gr+k 2 \f=Pf, (18) 



k d r d (m' + l) 2 (m' + 1) ,., . 

— — — — — a 4- h 'f — 

r dr k dr r 1r 

or the HEm'j, with 



u _ // cos m '0 
n ° l " 1/ sin m'0 



and f tf* <"'* \ (19) 

\— f cos m 0/ 



where / satisfies 

d im' - U 2 lm' - -\\ 

g r + k 2 \f=(3 2 f. (20) 



k d r d (m'-l) 2 . (m'-l) . ,,. _ y . 



/ dr k dr r 2r 

With m = m' ± 1, eqs. (18) and (20) can be combined as 
'k d r d 



rdrktr 7 * Tr 8 > + ' ^ '" '« 

The modes are indexed by the angular parameter m. For m = 
there are two polarizations for each ##!„ mode. This group includes 
the fundamental mode HEu, which propagates in single-mode fibers. 
For m = 1 there are two polarizations for the HEtn modes, also the 
TM modes and the TE modes. For m > 1 there are two polarizations 
for the HE m +i t „ modes and two for the EH m - Xn modes. 

When gradient index terms are neglected, eq. (9) reduces to the 
scalar wave equation 

For m = each solution represents two polarizations, i.e., two modes; 
for m ¥" each solution represents two polarizations and two angular 
harmonics, m' = m ± 1, i.e., four modes. 

The scalar wave equation and its counterparts [eqs. (10) and (21)] 
assume the form of a time-independent Schroedinger equation in two 
dimensions. In particular, the scalar wave equation can be expressed 
as 

r Jr r Jr ~ 7" + ( ** " **7 = ^ ~ k « )f > (23) 

where 

k cX = 27rn cl /A. (24) 

The quantity (k 2 — k 2 \) plays the role of a potential that is in the 
outer cladding and beyond; (p 2 — k 2 \) is an eigenvalue. 

Results on Schroedinger's operators 23 imply that the propagation 

OPTICAL FIBER PROPAGATION MODES 2669 



modes correspond on a one-to-one basis to the positive eigenvalues 
and the number of such modes is finite (see Ref. 23, p. 366). This fact 
implies that increasingly accurate estimates of the propagation modes 
can be obtained by ever finer discretization of the equations. By 
contrast, finer discretizations introduce ever more negative eigenval- 
ues corresponding to the continuum of radiation or unbound modes. 

2.2 Finite element reduction 

The Finite Element Method (FEM) can solve to desired accuracy 
the differential equations just derived. The present discussion special- 
izes to the scalar wave equation [eq. (23)]; modifications needed for 
the equations containing gradient index terms are indicated in the 
appendix. 

The solution function fir) is approximated by a piecewise linear 
function. This can be expressed in terms of interpolation functions 
(shown in Fig. 2) as 

L+l 

fir) = 1 f{rt)N,(r), (25) 

/=o 

where n = lb for / = 0, 1, • • • , L + 1 are evenly spaced sample points. 

The first and last terms of the series are affected by end conditions 
on f{r). At the center (r = 0) of the fiber fir) and its radial derivative 
must be bounded and well defined; so /(0) = when m > and df/dr 
(0) = or equivalents /(0) = fid) to first order in 5 when m = 0. 

At the other end, R = L5 represents a truncation radius in the 
cladding, and R + 5 represents the truncated part in a way that will 
now be explained. Assuming that the cladding extends indefinitely, 
the solution there is 

f(r) = aKMr), (26) 

where K m denotes the mth order modified Bessel function of the 
second kind, 24 a is an unknown coefficient, 

v = (/? 2 - &) 1/2 , (27) 

and m denotes the azimuthal mode number used in the scalar wave 



«//(/•)—. 




k--« — I 

Fig. 2 — Linear interpolation functions: hat functions. 
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equation [see eq. (23)]. It follows that the end condition for / at 
truncation is 

f'/f=vKL(vR)/K m ( V R). (28) 

A Taylor's expansion of / to first order about r = R yields 

f(R + 5) = f(R) + 8f'(R) 

= f(R)(l + 5rjK^r,R)/K m (vR))- (29) 

Most modes of a multimode fiber are unaffected by the end condition 
in eq. (29) because their waveforms are negligibly small beyond the 
core. But the waveform of a single-mode fiber can extend beyond the 
core and then the end condition needs to be enforced. 

Sample values of f{r), the coefficients in eq. (25), are estimated by 
the Galerkin weighted residual method. 22 Although a piecewise linear 
approximation cannot satisfy the scalar wave equation, it can satisfy 
weighted averages of the equation. In the Galerkin technique the 
weightings are chosen as the basis functions Nj(r) j — 1, • • • , L. In 
terms of inner products, defined for any functions p{r) and q(r) as 

J"R+5 
p(r)q(r)rdr, (30) 





- 3.^ -irfft. *& + <*• -«>/.*, 



Galerkin's technique yields 

-0* 1 -*&)(/, A» (31) 

for j = 1, • • • , L. The first term comes from an integration by parts. 

Substituting the piecewise linear approximation of / [from eq. (25)] 
into the Galerkin equations gives L equations for L + 2 sample values, 
but the end conditions eliminate two of these values. The result is L 
equations in L unknowns, expressed as the matrix equation, 

(A - m 2 B + C)f = 5 2 (0 2 - k\ x )Dt (32) 

The column vector f denotes the sample values, 

f = [/(ri), • • • , f{r L )t (33) 

The Lx L matrix A has jl element 

(dN t dN\ 

for all j and / except, to accommodate the end conditions, aio must be 
added to a u when m = and o^l+i/l+i/Zl must be added to a LL for all 
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m. The Lx L matrices B, C, and D have jl elements 

1 



c„ = ((k - k 2 cl )N h Nj)8 2 , (35) 

and 

djt = (N,, N,). 

The latter three matrices are evaluated by "lumping the masses," 
meaning that the integrals are evaluated numerically and the integra- 
tion points are chosen as the sample points. 25 The trapezoidal rule 
then yields for the off-diagonal elements and for the main diagonal 
gives 

b tt = 8 2 /l c u = l[k 2 (r,) - kl]b A da = Id 2 (36) 

for / = 1, • • • , L. The matrix A can be evaluated exactly. It is symmetric 
and tridiagonal (i.e., a jt = if | / - j \ > 1) and has 

au = -215 2 ay +1 - (/ + V 2 )5 2 (37) 

for / = 1, • • • , L. When m = 0, a n is -1.55 2 . 

Eq. (32) converts to a standard symmetric eigenvalue problem by 
multiplying both sides by the diagonal matrix D~ 1/2 and putting 

g = D 1/2 1. (38) 

The result is 

Tg - D~ l/2 (A - m 2 B + C)ZT 1/2 g = & 2 {(3 2 - k 2 cl )g, (39) 

standard form relative to the vector g. 
The key matrix, T, is symmetric and tridiagonal. For m * 

t„ = -2-y + 8 2 [k 2 (r,) - k 2 ci ] 



1 

ti+ij — tij+i — 



i + iV 2 ( i v 



(40) 



/ / \l + I, 

for / = 1, • • • , L - 1; for m = 

t n = -% + 5 2 [k 2 (r,) - /2c 2 .]; (41) 



for all m 



tu. = "2 - T2 + 5 2 [kHr L ) ~ fee 2 .] + t LtL+lgL+l /g L . (42) 

Li 



The end condition in eq. (29) combined with eq. (38) gives 
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gL + i/gL = —j— [1 + 8 v K^( v R)/K m ( v R)]. 



(43) 



Equation (39) represents the FEM reduction of the scalar wave 
equation. The other differential equations, as indicated in the appen- 
dix, reduce to the same form, with T symmetric and tridiagonal. The 
sample distance 8 is discussed in Section III. 

2.3 Calculation of the propagation modes 

The propagation modes are calculated by solving for the positive 
eigenvalues (//) of T. The associated propagation constants are 

(8 = W5 2 + kltf' 2 (44) 

with effective index of refraction 

n e = 0/(2ttA) = 0(c/o>). (45) 

When the end condition is enforced, the T matrix depends on fi in its 
(L, L) entry, but a simple iteration procedure yields the proper /?. The 
associated waveforms are given by f or g. 

The modal delays per unit length (r g ) are the reciprocal of the group 
velocities, 



dp 



(46) 



An efficient calculation of these uses the formula, 



d/3 






= - 



dk c \ 




1 


(dT \ 


' 


da) 2 


+ 


8 2 


WV 


g 



(47) 



where g is assumed to be normalized (i.e., g-g = 1). This follows 
standard procedure for taking the derivative of an eigenvalue with 
respect to a parameter. 26 Equations (40), (41), and (42) for T imply 
that dT/du 2 is a diagonal matrix. Using the equivalence between co 2 d/ 
die 2 and -\2d/d\ 2 , eq. (47) becomes 



L 
1=1 



*h> -x'gw 



gf/cn e . 



(48) 



To form the T matrix, the index profile n(r) must be available for 
any excitation wavelength (X). Sellmeier expansions 27 of the form 



n 2 = 1 + I Ai/[1 - (/,/X) 2 ] 



(49) 



have been fitted to measured values of refractive index over the range 
of wavelengths 0.8 ^m to 1.5 /*m for bulk samples of pure Si0 2 , 13.5- 
percent Ge doped Si0 2 , and 1 percent F doped Si0 2 (denoted here as 
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a, b, and c, respectively). Also, results in Ref. 28 indicate that the 
index is approximately linear with concentration. Therefore, the index 
profile of an optical fiber having a graded Ge dopant is taken as 

n(r, X) = a(X) + C„(r)[6(X) - a(X)], (50) 

where C n (r) denotes the concentration profile for Ge, relative to 13.5 
percent. For dual dopants, Ge and F, the index profile is 

n(r, X) = a(X) + C n (r)[b(\) - a(X)] + C n (r)[c(\) - a(X)], (51) 

where C n (r) denotes the concentration profile for F, relative to 1 
percent. 

The concentration profiles in these equations may be specified or 
may be deduced from the index profile at a reference wavelength. Once 
the concentration profiles are available, the index profile and its X 2 
derivative can be determined from eq. (50) or eq. (51) for any wave- 
length. 

The T matrix can be expressed in dimensionless form by replacing 
the sample spacing by 

5 = R^/U, (52) 

where R x denotes the core radius and Li the number of samples in the 
core. This gives 

5 2 [k 2 (r) - kl] = (2irR 1 /\) 2 (n 2 (r) - n 2 ,)/L 2 , (53) 

and from eq. (50), (with a = n ci ), 

n 2 (r) - n c 2 , = C n {n\ - nl) + (C 2 - C„)(m - n cl ) 2 , (54) 

where now C n (r) represents the Ge concentration normalized with 
respect to n lf the index at the center. The latter term in eq. (54) can 
usually be neglected because rii ~ n c i and often, C n (r) ~ 1 for most r. 
Then the T matrix can be expressed in terms of the usual V number, 

V = (27ri?x/X)(n? - nl) 112 = R x {k\ - fc c 2 ,) 1/2 (55) 

and the effective V number, 

V e = (2-KR 1 /X)(n 2 - n 2 ci ) 1/2 = R 1 W 2 - kl) 1 ' 2 (56) 

for the (L, L) element. The FEM reduction becomes 

T(V, V e )g = (Ve/Lrfg, (57) 

where now an iteration on V e is required. A similar expression, involv- 
ing three V numbers, holds for dual dopants. 
The effective index (n e ) is obtained from V, in eq. (56). The delay 
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IK 



T g 


_ w dp 2 _ 

" p du> 2 


1 dV^dV 2 

ff? dV 2 du> 2 + 


da> 2 




1 1 - 
cn e 


2 dfi& V e 





(rc? - n 2 ,) - X 2 —2 (rcf - n 2 i) 



(58) 



Matrix T is always symmetric and tridiagonal. The routine BISECT 
in the EISPACK 1 library computes the positive eigenvalues for such 
matrices and TINVIT, also in EISPACK, computes the associated 
eigenvectors g. These routines operate efficiently and reliably. 



2.4 Leaky modes 

When the index profile drops below n c \ over part of its range, then 
leaky modes may exist. These modes have n e < n c \ and complex prop- 
agation constants. The corresponding attenuation accounts for radia- 
tion loss as the waveform spreads out radially. As is well known, leaky 
modes also can arise when m # from the negative term, — (m 2 /r 2 ), 
in the propagation equation. 

Leakage loss is calculated by truncating the waveform at the begin- 
ning of the outer cladding and enforcing the proper end condition. In 
terms of the complex propagation constant y, eq. (27) for tj changes 
to 

i? = (-T 2 " *c 2 .) 1/2 , (59) 

and the end condition [in eq. (43)] becomes complex. Now, the T 
matrix has a real (T r ) and an imaginary (T/) part, but matrix T/ 
consists of all zeroes except the (L, L) diagonal element (call it x). 
First-order perturbation theory estimates an eigenvalue of T as 

n = n r + ixgl = 5V, (60) 

where Hr is an eigenvalue of T r and g L is the Lth component of the 
associated normalized eigenvector g. Again, an iteration is required 
because T depends on y. 



III. EXAMPLES OF MODAL CALCULATIONS 

In this section modal calculations based on the results of Section II 
are illustrated in 10 examples. The next section summarizes and 
evaluates the results. 
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3. 1 Infinite parabolic profile 

The parabolic profile that extends without radial limit has index 

n(r) = n x [l - 2A(r/J?) 2 ] 1/2 (61) 

for all r, where R denotes the nominal radius of the fiber core. Under 
the scalar wave approximation, its modal waveforms are Laguerre 
Gaussian functions, its propagation constants are 

ft«n = [kWi ~ 2Qk cini (2Ay /2 /R] 1/2 , (62) 

and its modal delays (neglecting material dispersion) are 

r m = (& m + kWi/(3» m )/2u>, (63) 

where the mode number is 

Q = 2 M + m + 1 m, m = 0, 1, ■ • - , (64) 

with m the angular harmonic and /i the radial harmonic. Each Q 
represents a group of modes that have the same propagation constant 
and delay (see Ref. 9). 

Modal delays were calculated assuming parameter values A = 0.013, 
R = 25 jttm, and X = 1.3 [im. Convergence with respect to truncation 
radius (TR) was complete (within 0.1 ps/km) for each mode for TR = 
1.512. Convergence with respect to the number (Li) of sample points 
in the core was within 1 ps/km in the rms delay for L x = 400, as 
indicated in Table I. An accuracy of 1 ps/km determines the bandwidth 
within 5 percent for a 10-GHz X km fiber and 0.5 percent for a 1-GHz 
X km fiber. The percent errors of the computed delays (with L x = 500) 
are shown in Fig. 3. 

The most accurately computed mode within a mode group had radial 
number n = 0; their waveforms have no zero-crossings and are most 
easily approximated by piecewise linear functions. The least accurately 
computed mode (which was off by 5.5 ps/km) had the largest fi. As 
another measure of accuracy, the spread in the computed delays for 
each Q is also given in Table I. 

To test the single-mode case, the parameters were changed to A = 

Table I — Convergence results for 

infinite parabolic profile in units of 

ps/km* 



In 


TRMS 


Spread (1) 


Spread (2) 


100 


144.0 


136.1 


96.5 


200 


129.4 


33.9 


23.7 


300 


127.1 


15.1 


10.5 


400 


126.3 


8.5 


5.9 


500 


125.9 


5.5 


3.8 



* (A = 0.013, R = 25 /tm, X = 1.3 |im) 
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Fig. 3 — Computation errors in modal delays for infinite parabolic profile. 



0.005 and R = 5 /xm, but again X = 1.3 fim. The calculated effective 
index of the HE n mode was accurate to seven decimal places and the 
delay to 0.1 ps/km when TR was 2.5i? and L\ was 300. This precision 
translates to nm accuracy in the zero-dispersion wavelength. The 
beam radius (BR) defined by the condition 



f 2 (BR)=f 2 (0)/e 



(65) 



was accurate to four decimal places, as indicated in Table II. The 
number of samples in the core can be less in the single-mode case 
because the HE n waveform does not oscillate; the truncation radius 
needs to be greater because the waveform extends farther into the 
cladding. 

3.2 Parabolic fiber 

Power-law fibers have the index profile 



n(r) ~ |n,(l - 2A) 1/2 



r<R 
r>R. 



(66) 
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Assuming parameter values of A = 0.013, R = 25 /tm, X = 1.3 jim, and 
a = 2 (the parabolic case), effective index (n e ) and delay {r g ) were 
computed for each propagation mode under the scalar wave approxi- 
mation. The rms delay (t RM s) converged to 1.793 ns/km within 0.3 
percent and the rms delay with the two highest-order-mode groups 
deleted (t rms) converged to 113 ps/km within 0.9 percent when L x = 
300 and TR = 0.1R. 
Figure 4 shows n e and r g for each mode arranged by mode group. 

Table II— Modal quantities for Ht u mode of parabolic 
profile* 



U 



TR 



n e 



(/is/km) 



300 3R/2 

300 2R 

400 2fl 

300 5ft/2 

Actual Values 



1.4531412 
1.4531608 
1.4531608 
1.4531609 
1.4531609 



4.8581500 
4.8577058 
4.8577059 
4.8577024 
4.8577025 



(A = 0.013, R = 25 Mm, X = 1.3 /im) 



BR (/xm) 



2.6561930 
2.6644067 
2.6644135 
2.6644397 
2.6643516 




4.898 



4.896 - 

> 

3- 4.894 

>- 

ui 

<=> 4.892 
ill 

□ 
o 
S 

4.8901- 



xxxxxxxxxxHJ«>< 

8 

X 
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X 
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MODE NUMBER 

Fig. 4 — Plots of n e and t 4 versus mode number for parabolic fiber. 
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The delays spread almost 10 ns/km for Q = 14 and 3 ns/km for Q = 
13, but less than 1 ns/km for the others. The spread in n e is negligible, 
reaching a maximum of 0.01 percent. 

Figure 5 shows n e and t s plotted jointly on an expanded delay scale 
with outliers excised. To compare with WKB results, n e and t b for the 
corresponding infinite parabola are indicated. The reference delays 
exceed all in their mode group: by as little as 0.2 ps/km for each of 
the first nine mode groups, but by more than 34 ps/km and 86 ps/km 
for Q = 13 and 14, respectively. 

3.3 Contribution of gradient index terms 

Gradient index terms cause modal delays to spread still further. 
Figure 6 shows the delay differences caused by these terms for each 
mode of the parabolic fiber. Of the 112 modal delays only 14 differed 
by more than 20 ps/km from the corresponding scalar wave values. 
The TM mode nearest cutoff had the largest difference (33 ps/km). 
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Fig. 5— Plots of T g versus n e for cladded and infinite parabolic profiles. 
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Fig. 6— Absolute differences in modal delays due to gradient index terms for parabolic 
fiber. 



The change in t RM s was 5.9 ps/km and in t'rms it was less than 1.5 
ps/km. The error is less than 1.5 percent in either case. In the examples 
that follow gradient index terms are neglected. 



3.4 Bandwidth of the parabolic fiber vs. X 

Material dispersion causes the bandwidth of power-law fibers to 
depend sharply on the excitation wavelength (X). Here and in the 
following examples, the profile is assumed known at the reference 
wavelength of X = 0.8 /im; the concentration profile and the index 
profile at other wavelengths are determined, as discussed in the 
previous section. 

Bandwidth is defined as the half -power frequency of the transfer 
function over some distance. When intramodal dispersion is neglected, 
the transfer function is 
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G(co) = £ a n exp(iW„), 



(67) 



where t„ denotes the delay and a 2 n the power of the nth mode. Defining 
the rms delay (t RM s) by 



trms — 
and the average delay by 



X a n (T n - r av ) 2 / X a n 



1/2 



T av — 2j &nJn / 2j ^n 



(68) 



(69) 



then to second order in the w(t„ — T av ), the bandwidth is 

BW= 1A27TTRMS), 



(70) 



in units of GHz x km for delays in ns/km. 

Bandwidth was calculated on this basis over a range of wavelengths 
for the parabolic fiber of example 2, assuming the tapered modal- 
power distribution shown in Fig. 7. Lower-order modes are weighted 
more heavily than the higher and the highest are omitted as the higher 
modes are increasingly vulnerable to microbending and cladding ab- 
sorption. 

Figure 8 shows the spectral plot of bandwidth. The peak value 
occurs at A =* 0.963 /*m, compared to 0.983 i*m reported in Ref. 13. 

The spectral plot is also shown for the case where material dispersion 
is neglected. The figure shows that the bandwidth is lower and essen- 
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Fig. 7 — A plot of the tapered modal-power distribution. 
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Fig. 8— Plots of bandwidth versus wavelength for parabolic fiber with and without 
material dispersion, using tapered modal-power distribution. 

tially constant except for small discontinuities associated with new 
modes cutting in. Other calculation indicates that bandwidth increases 
to a peak at 9.5 GHz X km for a « 1.98 compared to an optimum a « 
1.97 estimated by WKB analysis. 4 

3.5 Optimum a for power-law fibers 

Material dispersion causes the optimum a of power-law fibers to 
depend on excitation wavelength. Figure 9 shows bandwidth versus a 
when A = 0.013 and R = 25 urn for X = 0.82 fim and X = 1.32 ^m. The 
optimum a's are 2.07 for the former and 1.88 for the latter, compared 
to values of 2.081 and 1.884, respectively, reported in Ref. 13. The 
peak bandwidths are 5.48 GHz X km and 6.02 GHz X km, respectively. 

3.6 An alternate modal-power distribution 

The fractional power that propagates in the cladding can be used to 
derive an alternate modal-power distribution. Modes with more than 
0.1 percent of their power in the cladding will lose most of their power 
over 1 km for typical cladding losses of 1 dB/m. A realistic modal 
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Fig. 9 — Plots of bandwidth (in power-law fibers) versus a for X = .82 pm and X 
1.32 fim, using tapered modal-power distribution. 



distribution would neglect these modes and for simplicity could weight 
the others equally. 

Figure 10 shows bandwidth versus a for the two cases of example 5, 
assuming the new distribution. The optimum a is lower by 0.01 and 
the peak bandwidth is increased by 40 to 50 percent in either case. 
Table III indicates the modes neglected in the bandwidth calculation 
when X = 1.32 ^m and a = 1.88. 

3.7 Layer structure 

Layer structure in actual fibers perturbs the ideal profile. In this 
example power-law profiles are approximated by steps that span equal 
areas and match the ideal profile at the mid-area points of the layers 
(see Fig. 11). Bandwidth at X = 1.32 /xm (using the modal distribution 
in example 6) is shown versus a and the number of steps (NS) in Fig. 
12. The graphs indicate that the optimum a is essentially independent 
of NS, but the peak bandwidth and its sharpness decrease for smaller 
NS. 
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Fig. 10— Plots of bandwidth versus a for A = .82 nm and X = 1.32 /tm assuming equal 
excitation of modes having less than 0.1 percent power in cladding. 



3.8 Delay vs. X for a single-mode fiber 

Delay (per km) of the HE n mode was computed versus X for an 
actual single-mode fiber based on its measured index profile. The 
profile measurement was performed on the associated preform and is 
shown in Fig. 13. The depression of the inner cladding is due to F 
doping. 

As the calculation presumes radially symmetric profiles, the right 
and left sides were considered separately. About 750 sample points 
were used in either case. 

The radial scale of the fiber profile was assumed to depend linearly 
on that of the preform. The fiber core radius (R ) was estimated from 
the preform core radius (RP), fiber diameter (DF), and preform 
diameter (DP) as 



R = 1M(RP)(DF/DP), 



(71) 



where the 1 -percent addition estimates the Si0 2 loss in the outer 
cladding during the draw process. 
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Table III — Modes with less than 

0.1 percent of their power in 

cladding* 





<0.1 Percent 


>0.1 Percent 


Q 


m 


/* 


m 


M 


14 


1 

3 


6 

5 






13 



2 
4 
6 
8 
10 


6 
5 
4 
3 

2 
1 


12 





12 


1 
3 
5 
7 
9 


5 
4 
3 
2 

1 


11 





11 





5 


8 


1 




2 


4 


10 







4 


3 








6 


2 






10 


1 


4 


3 
5 


3 
2 








7 


1 








9 






* (A = 0.013, R = 25 urn, X = 1.32 /im, 
a = 1.88) 



Figure 14 shows computed differential delay versus X for the right 
and left profiles together with measurements. It was assumed that 
measurement matched calculation exactly at X = 1.32 nm, because 
only relative delays were measured. The zero-dispersion wavelength 
(at the delay minimum) was computed as X = 1.3127 /tm for the left 
profile and X„ = 1.3113 /im for the right, compared to measurement of 
X„= 1.3114 ^m. 

3.9 Leakage loss vs. X for a single-mode fiber 

Leakage loss (in dB/m) of the TE mode of the single-mode fiber in 
the preceding example was calculated for both sides of the profile as a 
function of X and is shown in Fig. 15. The difference between the left 
and right profiles becomes evident in this figure. A loss of 4 dB/m has 
been identified 29 with effective cutoff, and corresponds to cutoff wave- 
lengths, X c = 1.17 urn for the left profile and X c = 1.23 /xm for the right, 
compared to a measurement of 1.28 /xm. This discrepancy is discussed 
in the next section. 

The equivalent step approximation (where the index of the core and 
the inner cladding are area-weighted averages of the profile) gives 
X = 1.308 pm and Xc = 1.21 nm. The X c value straddles the previous 
calculated values, but the X value is somewhat less. 
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Fig. 11— A plot of a multilayered fit to an a profile (NS = 5, a = 1.90). 

3.10 Design curves for single-mode fibers 

The dimensionless formulation allows design curves to be generated 
efficiently. Figure 16 shows V e versus V for the HE n mode of three 
power-law profiles. The specific parameters (R, A, X) determine V and, 
hence, V e and dV e /dV. Delay versus X is computed according to eq. 
(58) and X is estimated at the minimum of the delay curve. Figure 17 
shows X versus R for two values of A for triangular profiles (a = 1). 

IV. SUMMARY AND CONCLUSIONS 

A method for computing modes of a circular optical fiber has been 
presented. The finite element method reduces Maxwell's equations to 
the standard eigenvalue problem, involving tridiagonal matrices. Rou- 
tines from EISPACK exploit the tridiagonal form to compute the 
eigenvalues and eigenvectors efficiently. From these the modal quan- 
tities are obtained. 

Using a piecewise linear approximation of the waveform is necessary 
to get tridiagonal form. Although piecewise quadratic and cubic ap- 
proximations of the waveform can lead to smaller matrices, tridiagonal 
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Fig. 12 — Plots of bandwidth versus a at X = 1.32 ^m for four values of NS. 

form is lost and the eigencalculations would be less efficient. Extension 
to elliptical and other nonradially symmetric fibers leads to similar 
difficulty in the eigen calculations. 

The method applies to any circular fiber and can account for 
gradient index terms to first order. As illustrated in the 10 examples 
of Section III, the calculations can provide information on the radial 
distribution of propagating power, on pulse dispersion (arising from 
material, intermodal, and intramodal dispersion), and on leakage loss. 

The accuracy of the method was tested for the infinite parabolic 
profile of the first example. After convergence was established, the 
modal quantities were compared with the actual values known through 
analysis. Agreement was excellent for both a multimode and a single- 
mode profile. 

The WKB and scalar wave approximations were tested for the 
parabolic fiber in the next two examples. The cladding of the fiber 
altered the delays of the higher-order modes substantially, a facet 
missed by WKB analysis. The gradient index terms contributed less 
than 1.5 percent to the rms delay; and so, for most purposes the scalar 
wave approximation suffices for this fiber. The validity of these 
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Fig. 13 — The preform profile for a single-mode optical fiber. 

approximations may change in other fibers; the value of the method 
is that it permits the test. 

Material dispersion was incorporated into the calculation in the 
remaining examples. The variation of refractive index with wavelength 
and dopant concentration was modeled by Sellmeier expansions based 
on measurements of bulk samples with certain dopant concentration. 27 
A linear dependence of index on concentration was assumed, but is 
not essential to the method. Irreversible thermal and stress effects 30 
might also be incorporated. 

Bandwidth was estimated in examples 4 and 5 for power-law fibers 
assuming a tapered modal-power distribution. With the higher-order 
mode groups deemphasized, results essentially agreed with WKB 
analysis. The usual sharp peak in bandwidth occurred in the spec- 
tral plot for the parabolic (a = 2) fiber in example 4 and in the 
a-dependence in example 5. Calculation of optimum a agreed reason- 
ably well in these cases with other numerical work and WKB calcu- 
lations. 

The bandwidth calculation was repeated in example 6 using a 
different modal-power distribution. Modes with more than 0.1 percent 
of their power in the cladding were neglected, as being too lossy to 
maintain power. The optimum as were essentially the same, but peak 
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Fig. 14— Measurement and calculations for left and right profiles of differential delay 
versus wavelength. 

bandwidths were somewhat higher. This example could be expanded 
to simulate differential mode attenuation, mode mixing, concatenation 
of dissimilar fibers, and other longitudinal variations. 

Perturbations of a profiles usually lower bandwidth. Example 7 
concerned the effect of layer structure on bandwidth and showed the 
expected decline with a smaller number of layers. Other perturbations 
such as ripple can also be treated. As a measure of its efficiency, the 
method used about two seconds per profile on the Cray-1 for this 
example. 

The last three examples concerned single-mode fibers. The first two 
of these involved an actual fiber whose profile was measured in the 
preform stage. Calculation of delay (per km) versus X matched meas- 
urement extremely well, but calculated leakage loss exceeded meas- 
urement, giving a 7-percent average underestimate of cutoff wave- 
length (A c ). Other calculations of leakage involving only slightly de- 
pressed inner claddings, where X c matches measurement to within 1 
percent, suggest that the discrepancy may involve material changes in 
the F-doped glass caused by the draw process (consistent with obser- 
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Fig. 17— Plots of computed zero-dispersion wavelength (AJ versus core radius (/?) 
for two values of A, for triangular profiles. 



vations in Ref. 31). Such changes need to be understood when pre- 
dicting transmission performance from a preform profile. Similar 
values for X and X c were obtained for the equivalent step profile (with 
depressed inner cladding), but the lower X indicates the effect of 
profile structure. 

The dimensionless formulation permits the greatest efficiency in 
getting design curves. Zero-dispersion wavelength (A ) for triangular 
profiles is calculated in example 10 showing the expected shift 32 to 
higher wavelengths. Spot size or cutoff wavelength can be obtained 
with similar efficiency. 

In summary, the calculation method is reliable, and relatively in- 
expensive. In the context of circular fibers it is comprehensive, capable 
of simulating diverse effects. 
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APPENDIX 

Gradient Index Terms 

Gradient index terms, involving the quantity q r = d/drlnk 2 , occur 
in eqs. (10) and (21). Although the index or its derivative may be 
discontinuous, they can be approximated by smooth functions, so k 2 
and g r can be taken as well defined, continuous functions of r. This 
appendix concerns the changes needed in the T matrix to accommo- 
date these terms. 

Equation (10) converts to symmetric form when h e = kf and both 
sides are divided by k. The result is 



Applying the Galerkin technique gives for the first term, 
A' ~ a'„ = - 



h i m >- 1 <"» 



L (f.fH<^> 

+ i(*f.*) + |Uf 



(73) 



in place of matrix A specified in eq. (34). The first quantity is O/i as 
before, the second adds to the C matrix of eq. (35), and the last two 
when evaluated by the trapezoidal rule contribute only to the first side 
diagonals. The term g r /r is handled in the same way as k 2 . Combining 
these contributions gives a new T matrix that is symmetric and 
tridiagonal. 

Likewise, eq. (21) converts to symmetric form when k 1/2 f replaces /. 
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The rest proceeds in the same way to give a symmetric, tridiagonal T 
matrix. 

Material dispersion is incorporated by expressing the index in terms 
of concentration functions as before. Radial derivatives needed for the 
gradient index terms themselves involve derivatives of the concentra- 
tion function. The co 2 derivatives of these terms, needed for the modal 
delays, are found by differentiating the coefficients of the concentra- 
tion functions. 
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